Many studies demonstrate the impact of ecological factors, such as geographic range, group size, and habitat breadth on the anatomical and behavioral adaptations of species. Group size is typically defined as the number of conspecifics with which an animal spends the majority of its time, and there is some evidence that these individuals form a cohesive unit (e.g., Jones et al., 2009). Habitat breadth has been variably defined as the number of geographic or climatic habitats used, dietary breadth, or even abundance. In particular, geographic range size and its utilization often correlate with abundance of a species (Brown, 1984, Lawton, 1993 and Pyron, 1999). Brown (1984) argued that this relationship occurs via the association between geographic range and resource abundance. Specifically, species with larger ranges typically have access to and utilize a wider variety of resources, and are consequently more likely to become widespread and abundant (Brown, 1984). In the present study, we follow Jones et al. (2009) in defining habitat breadth as the “number of habitat layers used by a species measured using any qualitative or quantitative time measure, for non-captive populations”. While we recognize that there are many possible definitions for the term habitat breadth, Jones and colleagues’ PanTHERIA database contains data on habitat breadth (using their definition) for close to 3000 species, so for data compilation purposes in our study, it was the logical definition to employ.

Species occupying a broad habitat tend to be prevalent and abundant (Brown, 1984, Gotelli and Graves, 1996 and Pyron, 1999). They also tend to be dietarily and behaviorally versatile, able to exploit a wide variety of local resources (Brown, 1984 and Brown, 1995). Species with a narrow habitat, on the other hand, tend to be habitat specialists, capable of effectively exploiting a narrow range of local resources (Brown, 1995 and Pyron, 1999). These specialists are more vulnerable to environmental fluctuations or changes in resource availability (Brown, 1984). Factors influencing resource exploitation may influence various gastrointestinal characters, and our data enable us to test this hypothesis.

By definition, species with larger groups live in closer proximity to conspecifics, consequently increasing their susceptibility to communicable diseases. Avoiding pathogens is a key ecological pressure influencing how individuals of a given species will distribute themselves in space (Brown, 1995). Therefore, it seems logical that species with larger groups may require anatomical and behavioral adaptations to living in closer proximity to one another than species with lower population density. Related variables include population group size: number of individuals that spend the majority of their time in a 24 hour cycle together; and group size (or social group size): similar to population groups size, except that there is some indication that these individuals form a social cohesive unit. If the immune functions of the appendix have played a substantial role in its evolution, we might expect it to be correlated with these characters.

Ecologists often focus on several key ecological variables when attempting to explain empirical patterns of anatomical, behavioral, and physiological diversity of species. These variables can offer valuable insights into the underlying mechanisms of ecological diversification among taxa, and help researchers explain the ecological “success” of species sharing niche space. Here, we focus on three of these variables – geographic range size, mean group size, and habitat breadth – given the influence that these factors have been found to have on socioecology and anatomical adaptation across mammals.

In order to determine the morphological variation and distribution of cecal appendices in mammals, data were compiled from the literature on presence/absence, size, and morphological configuration of the cecal appendix (Table 1 and Table 2; SOM 1–2). Two discrete characters represent appendix presence; the first, “observed” presence (character 2), represents scoring based on reports dealing with the species included in this study that specified whether or not the appendix is present. The second, “inferred” presence (character 3), scores more taxa (420 instead of 337) by taking into consideration statements in the literature that specify appendix distribution in larger, supraspecific taxa. As described above, we define a cecal appendix based exclusively on morphology, as a close-ended projection that is clearly differentiated from the cecum by a change in diameter.

In addition, published data were also collated on factors that have been suggested to co-vary with appendix presence, including gastrointestinal and dietary characters, as well as other ecological variables (SOM 1–2). In sum, the data included 14 discrete (Table 1, SOM 1–2) and 28 continuous characters (Table 2, SOM 1, 3). The anatomical data collected included characters of the cecum. In particular, the shape and size of the cecum, and cecal histological variables (including apex thickness and concentration of lymphoid tissue, often used in the definition of the appendix) were collated (Table 1 and Table 2, SOM 1). Our cecal shape characters included the state “paired ceca or colonic appendages”, in which we included instances of two equally-sized ceca, and instances of a second smaller “vestigial” cecum (e.g., Mitchell, 1905) and so-called “colonic appendages” (Stevens and Hume, 1995) or “distal paired caeca” (Bjornhag et al., 1994). As an outpocketing of the intestine, it has been suggested that the cecal appendix may be affected by other gastrointestinal characteristics. Inclusion of numerous gastrointestinal characters in the present study, such as size of the colon and gastrointestinal capacity, enabled us to identify possible patterns of co-variation among gastrointestinal traits (Table 1 and Table 2, SOM 1).

In addition to assessing co-variation among gastrointestinal traits, we chose to include characters that allowed us to empirically test the above-stated issues regarding: (1) a potential link between gastrointestinal characters (including appendix presence or size), and ecological variables relating to the geographic distribution and concentration of conspecifics in space that would affect the transmission of communicable diseases, such as group size and habitat breadth; (2) the possible relationship between diet and cecal appendix presence or size, as well as other gastrointestinal characters (e.g., Darwin, 1871). We also tested a third hypothesis, that the appendix is a functionless evolutionary vestige, by assessing whether it has appeared significantly more times than would be predicted by chance alone (based on a comparison between the number of gains and losses of the appendix). In order to test whether resource use and its corollaries, including degree of conspecific interaction, were correlated with appendix size and presence, we assessed: population density, geographic range, habitat breadth (as defined by Jones et al., 2009); as the number of layers used among the following: above ground dwelling, aquatic, fossorial and ground dwelling), mean social group size, home range, group range, and activity pattern (Table 1 and Table 2, SOM 1–2). The latter was tested because nocturnal species tend to be less gregarious than diurnal ones, and consequently might be expected to be exposed to lower levels of communicable pathogens, as well as differing in resource exploitation.

Dietary characters incorporated to address the dietary hypothesis for appendix evolution included: diet and gut adaptation, cellulose content of diet, type of fermentation (e.g., foregut, hindgut, etc.), relative dry matter intake, particle retention time, food quality, and dietary breadth (i.e., number of dietary categories eaten by each species over a period of time) (Table 1 and Table 2, SOM 1). It has recently been argued that dry matter gut contents may be a more reliable indicator of an animal's ability to process low quality food sources than is mean particle retention time [MRT] (Müller et al., 2013); however, due to the relative paucity of data on dry matter gut contents across mammals, we chose to include MRT instead. Ecological characters that could potentially affect diet and digestive strategies were also considered, including precipitation, temperature, latitude, and terrestriality (Table 2).

The correlation test that we used (see below) requires that the states be ordered according to clines. This was often straightforward, but for cecal morphology, this was less intuitive. We ordered the states along a gradient from an appendix-like, slender to a bulky cecum. Even though either end of this spectrum could have been scored “0”, we selected as the state “small, appendix-like cecum” as being “0” because it is probably primitive for mammals, as shown by the fact that it is present in all sampled monotremes and some marsupials, including Caenolestes, which is the sister-group of the other sampled marsupials, in our reference tree (ordered parsimony unambiguously supports this interpretation). This character was scored as inapplicable for taxa lacking a discrete, recognizable cecum, both for logical reasons (the shape of an unrecognizable structure cannot be known) and because optimization of the character, with states ordered (as is most appropriate for characters with clines) or not, suggests that absence of a cecum results from a loss, and that this may occur from most states of this character. Thus, there is no logical place to insert the state “absent” into the shape morphocline. But given uncertainty about the proper ordering scheme between some states of cecal morphology, we performed a sensitivity analysis by lumping states and re-testing for correlation with the appendix; thus, the number of states was reduced from six (excluding absence, scored as inapplicable) to five, four (two alternative schemes examined), and three (Table 1).

It has been suggested that age and body size at time of weaning impacts the quantity of maternal gut bacteria that are transferred from mother to infant (Bezirtzoglou et al., 2011, De Leoz et al., 2012 and Martín et al., 2003). Species that nurse their offspring for longer periods should gain a higher concentration of gut bacteria (Bezirtzoglou et al., 2011 and Langer, 2003). Thus, we also tested whether weaning age and weaning body mass were correlated with appendix presence.

Data were collected for 533 terminal mammalian taxa comprising all speciose mammalian families (SOM 1–2). Two more taxa are included in the database, but the dog (Canis familiaris), and pacas (Cuniculus) were excluded in the first case, because the dog is a direct descendant of the wolf, which is also included in our study, and in the second case, because too few gastrointestinal and ecological data were available.

The phylogeny that we used in our previous paper (Smith et al., 2013) was obtained from Bininda-Emonds et al., 2007 and Bininda-Emonds et al., 2008. Here, we updated it (for topology and branch lengths), whenever we could find more recent information (which is the case for the vast majority of taxa) using Meredith et al. (2011) for all taxa included in the latter except for a few cases in which this contradicted the established consensus and other recent studies (see below). For several clades, Meredith et al. (2011) offered insufficient resolution (too few taxa were included). Thus, we also used: Phillips et al. (2009) for the divergence date between Zaglossus and Tachyglossus, Voss and Jansa (2009: fig. 35) for Didelphidae (topology only; this happened to be coherent with Bininda-Emonds et al., 2007), Meredith et al. (2008) for Macropodiformes, Meredith et al. (2009: fig. 3) for divergence times within Acrobatidae and Diprotodontia in general, Meredith et al. (2010: fig. 5) for Pseudocheiridae, Krajewski et al. (2000: fig. 4) for divergence times within Dasyurus, Delsuc et al. (2012) for xenarthrans, Moraes-Barros et al. (2011) for the phylogeny within Bradypus, Castro et al. (2013) for Dasypus (topology only), Poux et al. (2008) for Tenrecidae, Smit et al. (2011) for Macroscelididae, Ohdachi et al. (2006: fig. 1) for Soricidae, Esteva et al. (2010) for Sorex, Dubey et al. (2008) for Crocidura, Tougard et al. (2001) for Rhinocerotidae, Agnarsson et al. (2011) and Teeling (2009) for chiropterans, Prevosti (2010) and Agnarsson et al. (2010) for carnivoran topologies between low-ranking taxa (divergence times between these were not always available; in these cases, they were kept as close to those originally in Bininda-Emonds et al. (2007) as possible, if this was compatible with Meredith et al. (2011), and if no more detailed data were available in other studies), Slater et al. (2010) for ursid phylogeny, Patou et al. (2009) for Herpestidae, Yonezawa et al. (2007) for mustelids, Koepfli et al. (2008) for mustelids, Johnson et al. (2006) for felids, Bagatharia et al. (2013) for divergence times within Panthera, Hassanin et al. (2012: table 1, UNI-HARD mean column for ages) for Cetartiodactyla, Bibi (2013) for Bovidae and closely related taxa, Perelman et al. (2011) for primates, with some additions from Masters et al. (2007: fig. 4) for galagonids, and Cortés-Ortiz et al. (2003) for Alouatta, Honeycutt (2009) for high-level topology and divergence times among rodents, Nunome et al. (2007) for Gliridae, Mercer and Roth (2003) for Sciuridae, Vilela et al. (2009) for hystricomorphs, Hafner et al. (2007) for heteromyids, Belfiore et al. (2008) for divergence times within Thomomys, Jansa et al. (1999) for Nesomyidae and Mystromys, Lebedev et al. (2012) for topology within Dipodidae, Salazar-Bravo et al. (2013: fig. 2) for the topology of sigmodontine cricetids, with divergence times from Parada et al. (2013), Ventura et al. (2013: fig. 2) for Chilomys instans, Bradley et al. (2007) for Peromyscus, Jansa et al. (2006) and Lecompte et al. (2008) for Muridae, Galewski et al. (2006) for Arvicolinae, Neumann et al. (2006) for Cricetinae, Robovský et al. (2008: fig. 5) for Microtus, Steppan et al. (2004) for muroid topologies, Rowe et al. (2008) for Murinae, Steppan et al. (2005: fig. 6) and Colangelo et al. (2007) for Gerbillinae, and Verneau et al. (1998) and Robins et al. (2008) for Rattus. For a few taxa that were not covered in these studies, divergence times were obtained from Kumar and Hedges (2011).

For a few taxa, sources conflicted. In such cases, we opted to retain topologies that reflected the established consensus, rather than new suggestions (Fig. S1). Thus, we placed Scandentia as sister-group of the clade that includes Cynocephalidae plus primates, as recently upheld by Lartillot and Delsuc (2012), rather than placing it as the sister-group of Glires, as recovered by Meredith et al. (2011). We deliberately did not use automated supertree construction methods, for three reasons. First, such methods are aimed at giving an automated consensus when there is significant overlap between sources and when there is no a priori reason to prefer one source tree over another, in case of conflict. This is not the case here as overlap was minimal between sources, and the few conflicts could be easily resolved by using the most recent paper, which typically had used much more data. Second, there is a problem with such methods, such as MRP (Matrix Representation Parsimony), which may yield clades that are not included in any of the source trees (Goloboff and Pol, 2002). Third, such methods typically do not incorporate branch lengths, which are critical for some of our analyses.

We investigated the presence of a phylogenetic signal in the data to determine if phylogeny-informed tests were required (Laurin, 2004). Given that several of our characters are continuous, we used a test that uses squared-change parsimony to compare the amount of change implied by each character on the reference tree to the amount of change on a population of randomized trees. Given that squared-change parsimony is sensitive to branch lengths, the simplest way to produce such a population of random trees is to randomly reshuffle terminal taxa on the reference tree, whose topology (except for the identity of terminal taxa) and branch lengths are thus kept constant (Laurin, 2004). This procedure can also be used for discrete characters (as was done here), although other methods to produce random trees could also be used for these. These tests were performed in Mesquite (Maddison and Maddison, 2014).

Given that most characters display a strong phylogenetic signal (see results), we used only phylogeny-informed tests. We first checked if the continuous data could be analyzed through phylogenetic independent contrasts (PIC), the most commonly used method to analyze continuous comparative data in a phylogenetic context (Felsenstein, 1985). However, the four artifact checks implemented in Mesquite (Maddison and Maddison, 2014) indicated that most of these data departed strongly from a Brownian motion evolutionary model on the reference tree (Fig. S1). These checks consist of verifying the relationships between absolute value of standardized contrasts and expected standard deviation (based on branch lengths), between contrasts and nodal value, between contrasts and nodal height, and between nodal value and nodal height. For an ultrametric tree (i.e. with tips or terminal taxa all at the same height above the base node), none of these relationships should be statistically significant if the contrasts are adequately standardized. Most characters exhibited strong standardization artifacts (Table S1), even after correcting for multiple tests as suggested by Canoville and Laurin (2010). To attempt to correct this problem, the data were log-transformed, but this only moderately improved standardization, with most characters still displaying very highly significant artifacts (Table S2). Therefore, PIC were not used to assess correlations between the continuous characters; instead, the non-parametric, phylogeny-informed sign test known as (phylogenetic) pairwise comparisons (Maddison, 2000 and Read and Nee, 1995) was used to assess correlations between characters, continuous or discrete. This test is implemented in Mesquite (Maddison and Maddison, 2014). It consists of contrasting pairs of terminal taxa. We used mostly the pairing algorithm that contrasts taxa differing in the state of the independent character. However, the Mesquite implementation cannot handle missing data or polymorphism in the independent character, or polytomies in the tree, so the master tree was pruned to retain only terminal taxa that were scored (and non-polymorphic, and not with the “inapplicable” state). When the subtree retained a polytomy, it was resolved arbitrarily by Mesquite (without respect with character-state distribution). Contrary to our previous study (Smith et al., 2013), we did not randomly resolve polytomies ten times because our updated tree is much more resolved, with only a couple of small polytomies involving three taxa, so their impact on the analyses is negligible. In several cases, more than one pairing of terminal taxa is possible; in these cases, we averaged the probabilities of the first ten pairings found and report these values. Given the number of characters included in this study, not all possible pairs of characters were checked; instead, we focused on the possible relationships between the presence or size of the appendix and the other gastrointestinal anatomical characters, as well as between gastrointestinal anatomy and ecology.

To better assess the significance of the potential correlation between appendix presence and cecal morphology, we checked the proportion of taxa with each cecal morphology that possess an appendix. We also produced 2 × 2 contingency tables showing the distribution of character combinations in taxa and performed Chi² and Fisher's exact tests on these using Graphpad Quickcalcs (http://www.graphpad.com/quickcalcs/), to complement the pairwise comparison tests. Scatterplots are provided to support some of the results of pairwise comparisons between continuous characters (Fig. 2).

We compared the number of appearances and losses of the cecal appendix using parsimony optimization on the reference timetree. This test enabled us to empirically test the hypothesis that the appendix has appeared more frequently and disappeared less frequently throughout mammalian evolution than would be expected by chance alone if this character were selectively neutral. Rejection of the null hypothesis of equal rates of gains and losses would suggest that the appendix likely serves some adaptive purpose. Our results could differ substantially from those of Smith et al. (2013) because that previous study relied on a reference tree that had several large polytomies, which create problems with character optimization (Maddison, 1989). To deal with the ambiguity of optimization of appendix gain and loss on our tree, and considering that Smith et al. (2013) found that gains were much more numerous than losses, we used the most conservative number of gains (minimal number) and losses (maximal number) suggested by parsimony to compute the probability that both events are equally probable, which should be a conservative approach. This probability was computed using GraphPad Software (http://www.graphpad.com/quickcalcs/binomial2/). We used a two-tailed test because our initial hypothesis that gains are more probable than losses is derived from examination of the optimization, which is itself based on our data. Thus, using a one-tailed test would have introduced a bias into the study.

To compare evolutionary rates of the cecal appendix and habitat breadth between clades, we divided the number of transitions (appearances or losses) of the appendix by the sampled phylogenetic diversity (sum of branch lengths) of the clades (Faith, 1992), as Smith et al. (2013) did. We used a binomial test to determine the probability of the null hypothesis under which the proportion of events in each clade simply reflects its sampled phylogenetic diversity, which reflects the hypothesis that the appendix is equally likely to be gained (or lost) in all clades.

Given that comparative studies often include multiple tests, corrections are required when applying the probability threshold to assess statistical significance of the results. We applied the False Discovery Rate (FDR below) procedure described by Benjamini and Hochberg (1995) because it retains more power than Bonferroni corrections, while retaining statistical validity (Curran-Everett, 2000). In this case, we applied the procedure simultaneously on all our tests, including those pertaining to the presence of a phylogenetic signal (but not the presence of statistical artifacts linked with inadequate phylogenetic contrast standardization, because most of these were strongly significant and could have biased the analysis), as well as to the correlations between characters. This allowed us to assess the significance of 170 tests (and associated probabilities) simultaneously and to eliminate about 8–9 expected false positives (170 times 0.05 probability of obtaining a significant result yields 8.5 expected spurious results; actual number of initially significant results that become non-significant after correction for multiple tests is reasonably close, at 5).

The false discovery rate (FDR) procedure indicates that the correct probability threshold is 0.017 (instead of the 0.05 threshold that we would have used if we had not accounted for multiple tests).

All discrete (Table 3) and most continuous characters (Table 4) display a strong phylogenetic signal. The only continuous character that does not display such a signal is the relative dry matter intake (9).

Among discrete characters, appendix presence (observed or inferred) is significantly correlated with the concentration of lymphoid tissue, cecal morphology, and with cecal apex thickness (Table 5).

The correlation with cecal morphology is strongest when the latter is coded into six to four states; it is non-significant with only three states. The relative frequency of the appendix in taxa with the various states of cecal morphology is generally congruent with our hypothesis about state orders, with a few discrepancies (Table 6). The only exceptions concern the tapering cecum, which is more frequently associated with an appendix than the spiral cecum (though we had trouble ordering these two a priori based on their morphology) and the rounded cecum, which is more frequently associated with the appendix than the paired ceca/colonic appendages and cylindrical ceca (though again, ordering it compared with the paired cecum/colonic appendages state was especially problematic; Fig. 1). Yet, despite these mismatches, the relationship is significant, even after accounting for corrections for multiple tests, for more than half of the tested coding schemes, which suggests that a genuine relationship exists between cecum shape and appendix. This result is further confirmed by a classical, non-phylogenetic Chi² test (Table 6), whose strongly significant result (P < 0.0001) reflects the fact that the proportions of taxa with an appendix differ drastically between cecal morphologies.

The relationship between appendix presence and concentration of lymphoid tissue is supported by an examination of the distribution of the states, and by non-phylogenetic statistical tests (e.g., Fisher's exact test) on these (Table 7). Similarly, the relationship between appendix presence and cecal apex thickness is also supported by classical statistical tests (Table 8).

Tests between discrete (inferred appendix presence) and continuous characters yielded only one marginally significant result (with adult body length), which is no longer significant after corrections for multiple tests using FDR (Table 9).

Tests between continuous characters, namely length of appendix, cecum and colon (all divided by the cubic root of body mass to adjust for body size effects) and the other gastrointestinal and ecological characters, yielded several significant correlations, four of which remain significant after corrections through FDR (Table 10). Of the ecological variables of particular interest to the research questions, there is a positive correlation between cecum length and group mean size (P = 0.0013; Fig. 2), a negative relationship between cecum length and habitat breadth (P = 0.0073), and a significant inverse correlation between colon length and weaning age (P = 0.0101; Table 10). There are no significant correlations between appendix length and home range, dietary breadth, or any of the other dietary or ecological variables (Table 10). Among the gastrointestinal variables, unsurprisingly, we found a positive relationship between cecum and colon length (P = 0.0003).

Parsimony suggests a minimum of 29 gains and a maximum of 12 losses of the cecal appendix. Note that not a single loss event is unambiguously supported by our tree, so that there might actually be 41 gains and 0 losses. However, even the most conservative estimate of 29 gains and 12 losses, if we hypothesize that both events are equally probable, the test yields a probability of 0.0115 (two-tailed test) of finding such asymmetrical results. This leads us to reject the null hypothesis, and we conclude that gains in the appendix are more probable than losses.

These data support a great heterogeneity in evolutionary rates of the cecal appendix between clades (Table 11). The probability that the actual evolutionary rates are equal in Laurasiatheria (which show no transitions) and Euarchontoglires (which show 36 transitions) is inferior to 0.00001. Even the less spectacular difference in evolutionary rates between Euarchontoglires and Metatheria is statistically very highly significant (P = 0.00094).